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ABSTRACT 

We constrain the properties of massive binaries by comparing radial velocity 
data on 114 early-type stars in the Cygnus OB2 Association with the expecta- 
tions of Monte Carlo models. Our comparisons test several popular prescriptions 
for massive binary parameters while highlighting the sensitivity of the best-fitting 
solutions to the adopted boundary conditions. We explore a range of true binary 
fraction, F, a range of power-law slopes, a, describing the distribution of com- 
panion masses between the limits qi ow and 1, and a range of power-law slopes, /3, 
describing the distribution of orbital separations between the limits r in and r out . 
We also consider distributions of secondary masses described by a Miller-Scalo 
type initial mass function (IMF) and by a two-component IMF that includes a 
substantial "twin" population with M 2 ~ Mi. Several seemingly disparate pre- 
scriptions for massive binary characteristics can be reconciled by adopting care- 
fully chosen values for F, r in , and r out . We show that binary fractions F < 0.7 
are less probable than F > 0.8 for reasonable choices of r« n and r out . Thus, the 
true binary fraction is high. For F =1.0 and a distribution of orbital separations 
near the canonical Opik's Law distribution (i.e., flat; (3 = 0), the power law slope 
of the mass ratio distribution is a =-0.6 - 0.0. For F ~ 0.8, a is somewhat 
larger, in the range -0.4 - 1.0. In any case, the secondary star mass function is 
inconsistent with a Miller-Scalo -like IMF unless the lower end is truncated below 
~2-4 M & . In other words, massive stars preferentially have massive companions. 
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The best fitting models are described by a Salpeter or Miller-Scalo IMF for 60% 
of secondary star masses with the other ~ 40% of secondaries having M 2 ~ M 1; 
i.e., "twins". These best-fitting model parameters simultaneously predict the 
fraction of type Ib/c supernovae to be 30-40% of all core-collapse supernovae, in 
agreement with recent observational estimates. 

Subject headings: techniques: radial velocities - - (stars:) binaries: general 

- (stars:) binaries: spectroscopic — stars: early-type — stars: supernovae - 
gamma rays: bursts — X-rays: binaries — 



Introduction 



1.1. Significance and Theoretical Parameterizations of Binaries 



Massive binary stars are invoked as the progeni tors of a wide va riety of astrophysical 
phenomena, from s hort and long gamma -ray bursts (IFryer et al.lll999l ). type Ib/c and blue 
type II supernovae (jPodsiadlowskil Il992l ) . to the fastest runaway O/B stars (jBlaauwlll96ll ) 
and the entire menagerie of binary systems with compact r emnants: X- r ay bin aries, millisec- 
ond pulsar systems, and double neutron star systems (see IFryer et al.l (119981 ) for a review). 
Whether binaries are necessary for the formation of all such objects is a matter of debate. 
For instance, although many authors propose that type Ib/c supernovae and short- and long- 
duration gamma-ra y bursts are produc ed in binaries, it is still argued that their progenitors 
may be single stars (IHirschi et al.ll2005l ). Observational studies can contribute to this debate 
by measuring the frequency, separations, and mass ratios of close massive binary systems. 

Although debates on the role of binaries continue, it is generally agreed that most 
of the progenitors of X-ray binaries and double neutron star binaries in the Galactic disk 
must arise from massive binary systems. Population synthesis models of these systems, 
however, adopt very different properties for the progenitors' stellar components. These 
studies begin with four assumptions regarding the binary characteristics: the distribution of 
primary masses, the distribution of secondary masses, the distribution of orbital separations, 
and the distribution of eccentricities. Population synthesis studies all choose some analytic 
form for these distributions, generally basing them upon observational constraints on the 
initial and final binary system parameters. 

One significant source of disagreement in these initial distributions is the choice of the 
secondary star mass function. Typicall y, the primary star mass d istribution is chosen using 



a field initial mass function (IMF) (e.g.. lScalo 



19861 : ISalpeterlll955l ) . The secondary mass can 
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be chosen through the same IMF or, more commonly, from a distribution of the secondary- 
to-primary mass ratios: f(q), where q = M 2 /Mi. In this latter formalism, the distribution 
is chosen to be proportional to the mass ratio to some power: f(q) oc q a between the limits 
qiow < q < 1-0. Note that the power-law slope, a, is not equivalent to the power-law slope, 7, 
of the secondary star mass distribution, comm o nly ex pressed /(M2) oc . The relationship 
between f(q) and /(M 2 ) is discussed by iToutl (11991I ). and it is, in general, non-trivial. 



Population studies disagree on the choice of a, primarily because the choice of a strongly 
effects the formation rate of key astrophysical objects. High values maximize the formation 
rate of double neutron star binaries whereas low values maximize the formation rate of low- 
mass X-ray binaries. Genera lly, population synthesis stu dies of low-mass X-ray binaries have 
assumed a values near -2.7 (IKalogera &: Webbink!ll998l ). whereas studies of dou ble neutron 
star binaries tend to prefer a flat distribution with a = (IBelczynski et al.ll2002l ). Resolving 
this ambiguity is important because the choice of this value can change forma tion rates of, 
for example, neutron star binaries by nearly a factor of 100 (IFryer et al.l 119991 )! This is the 
difference between advanced LIGO detecting 10 binary inspirals per year versus 1 inspiral 
in 10 years. 

The second critical parameter is the distribution of orbital separations, f(r). This 
distribution is often taken to be a power law between some inner and outer limits, ri n 
and r out , such that /(logr) oc (logr) /3 where (3 > indicates a preference for large orbital 
separations and (3 < indicates a preference for small separations. 



1.2. Observational Constraints on Mass Ratios 



Observational studies of binaries have reac hed widely varying, even dis parate, conclu- 
sions r egarding these key parameters (reviewed in lAbtl (119831 ) : iLarsonl (120011 )). iGarmany et al. 
( 119801 ) found a = — 1 for a sample of 67 massive spectroscopic binaries, suggesting that 
the mass ratio distribution is ei ther flat or peaked toward massive companions. Among B2- 
B5 primaries , lAbt et al.l (jl990l ) report a secondary star mass distribution consistent with 
the ISalpeterl (jl955|)_ IMF, corresponding to a ~ —2.3. In a spectroscopic study of solar 
type field stars, (jDuquennoy fc Mayor! Il99lh conclude a < — 2 for systems with q > 0.1, 



consistent with secondaries drawn from a field star IMF. Adaptive o ptics imaging surveys 



for companions of B and A stars in the Scorpius O B2 Association ( jKouwenhoven 1 12006 



Kouwenhoven et al.ll2005l ; IShatsky fc Tokovininl 120021). however , find the mass radio distri- 

0) 



Fisher et al.l (120051 ) report a flat 



(a ~ 



bution to be well-characterized by a— -OA — -0.5 
distribution of mass ratios among single-lined spectroscopic binaries (SBls) in the solar 
neighborhood, while double-lined spectroscopic binaries (SB2s) show a peak near q ~ 1. 
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Arguments abound over whether this q — 1 peak is real or a selection effect arising from the 
fact that SB 2s are only identified as such i f the mass and luminosity ratios are near unity 
rtLucv fc Riccolll979l : [ifalbwachs et al.ll2003h 



Several investigations have concluded that there are two populations of secondaries 
leading to a bimodal dist ribution of mass ra t ios. In summarizing dat a from the Index Catalog 
of Visual Binary Stars (jJeffers et al.l 119631 ) IHogeveenl (I1990I. Il992l) and the 6th Catalog of 
the Orbital Elements of Spectroscopic Binary Stars (IBatten et al.lll989l ) found that a < —2 
for systems with mass ratios q > q where q ~ 0.3. This is roug hly equivalent t o secondary 



1955 



Scalo 1986 



masses being draw n from the field star initial mass function (iSalpeter 
Kroupa et al.lll99ll ). For extreme mass ratio systems, q < 0.3, IHogeveenl (119921 ) finds a 
flat distribution (i.e., a ~ 0). Furthermor e, the turnover point, q , may vary inversely 



with primary mass, but the IHogeveenl (119921 ) sample did not inclu de O stars. Using radia l 
velocity data of solar-type F7-K stars with periods P < 10 yrs, lHalbwachs et al.l (120031 ) 
find a broad peak in the mass ratio distribution near q = 0.5 and a sharp peak near q = 



0.8. In a spectroscopic study of low-mass stars, iGoldberg fc Mazehl (119941 ) report a peak 
in the mass ratio distribution near q = 0.2 and a sma ller peak near q = 0.8. From a 



spect r oscopic survey of eclipsing OB binaries in the SMC (IHarries et al.ll2003l ; iHilditch et al. 



20051 ) iPinsonneault fc Staneki (120061 ) conclude that massive primaries have two populations 
of companions: a "twin" population with M2 > 0.95Mi comprising 45% of binaries, and 
a population drawn from a flat mass distribution (f(q) = constant corresponding to a=0) 
comprising 55%. This report echos earlier findings among low-mass stars that short-period 
(P < 40 d) binaries have a "twin" com p onent compr ising 10% - 50 % of the total pop ulation 
(ILucy fc Riccol Il979l ; iTokovininl |2000| ) . iLucyl (120061 ) reassessed the IHogeveenl (Il992l) study 



using data from the 9th The Ninth Catalogue of Spectroscopic Binary Orbits (IPourbaix et al. 
20061 ) and concluded that the data support an excess of q ~ 1 "twin" systems. 



1.3. Observational Constraints on Orbital Separations 



Duquennoy fc Mayor! (119911 ) use a radial velocity survey of G dwarf primaries to con- 
clude that the distribution of orbital periods for nearby binaries is broad and approximately 
Gaussian in logP having logP = 4.8 and o"i og p = 2.3 with P in days. Therefore, the dis- 
tribution of orbital separations should be approximately log-normal as well, with a mean of 
~ 25 AU . This result contrasts with others that indicate the distribution of orbital separa- 
tions isjlat i n log space, f(l o gr) oc (log r)°, corresponding to f ir} oc (r)" 1 (i.e., Opik's law; 
6pikill924h . iKouwenhovenl (hood ) and IShatsky fc Tokovininl (120021 ) measured separations 
r > 30 AU for visual binaries in Sco OB2 and find a flat distribution consistent with Opik's 



law. ILepine fc Bongiornd (120071 ) report a flat distribution of separations among Hipparchos 
common proper motion stars in the solar neighborhood. 



1.4. Observational Constraints on Binary Fraction 



Estimates for the binary fraction, F, among massive stars vary widely, and such esti- 
mates are always lower limits given that very close, very distant, and very low-mass compan- 
ions are difficult to detect. Here, we use F to mean the fraction of all primary stars having 
exactly one secondary. We do not explicitly consider triple and higher order systems, even 
though iGarcia fc Mermilliodl (120011) r e port t hat a signif i cant f r action of mass i ve clo se bina- 
ries may have multiplicity > 2. kbtl Jl983h : kbt et al.l f|l99Clh : IPoveda et all (jl982h report 
binary fractions of at least F = 0.7, and suggest that the binary fraction may be clo se to 
100% after corrections for incompleteness. Among O stars in clusters and associations. iGies 
(119871 ) report binary fractions of F =26% - 55%. iMason et al.l (1 19981 ) find binary fractions 



of at least F =59% - 75% for O stars in clusters and associations. Garcia fc Mermilliod 



(120011 ) tabulate multiplicity statistics for O stars in a variety of cluster environments, and 
they conclude that the binary fraction varies between F =14% and 80%, with denser clusters 
having lower binary fractions. 



1.5. Scope of this Paper 



Despite a wealth of observational data on parameters of binary systems, there is little ap- 
parent agreement, leaving reasonable doubt as to whether the results are reconcilable. Much 
of the disagreement can be explained by observational and sample selection biases. Both 
spectroscopic radial velocity surveys and direct imaging surveys for companions of massive 
primaries become incomplete at low mass ratios below q ~ 0.25. Secondaries in close, extreme 
mass ratio systems with massive primaries are presently undetectable with either technique 
and are likely to remain so in the near future. Thus, estimates of F, f(r), and f(q) at low q 
must be made by extrapolating some functional form of the distribut ion based on da t a from 
systems w i th lar g er q and/or large r. Radial v el ocity surveys such aslGarmanv et al.l (|19 8Q N ) ; 



Abt et al 



( 119941 ): 



Harries et al. 



Fisher et al. 



preferential 



Jl990h: IPuquennov & Mavorl dl99lh Jgalbwachs 6^ (|2003h : iGoldberg & Mazeh 



companions (a < few x 10 AU) while imaging su r veys such as 



y detect close, short - perio d 



Shatskv fc To 



covininl (j2002l) : 



Kouwenhoven et al.l (120051 ) ; iKouwenhoven I (120061 ) ; ILepine fc Bongiornd (120071 ) preferentially 
detect more distant, long-period companio ns. If the fundamental parameters of binaries are 



different for close versus distant binaries (I Abt et al.l Il990l ; iMazeh et al.l Il992l ; IMason et al. 
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19981 ). then the variety of reported results might be partially under stood in ter ms of these 
observational selection effects (see Fig 1 of iMason et al.l (119981 ); also iToutl (ll99ll )). Further- 
more, many studies do not explicitly describe the limits, r in , r out over which the results are 
valid, complicating comparisons between surveys. 

The statistics on companions to massive stars ar e especially uncertain as there have 
been few studies of massive binaries since the survey of iGarmany et al.l (Il980l ). For massive 
primary stars, the small number of broad spectral feature s, coupled with r otatio nal broaden- 
ing, make radial velocity measurements difficult, and the lGarmany et al.l (119801 ) sample may 
have missed binaries below velocities of 50kms _1 . More than half of their sample consisted 
of double-lined spectroscopic binaries, thereby introducing a bias in favor of nearly equal 
mass components with large relative radial velocities. About half of the sample consisted 
of evolved O stars, selected on a magnitude-limited basis from throughout the sky. Given 
the selection biases an d lack of a common origin and age for this sample, the conclusions of 
Garmany et al.l (119801 ) may not apply to the general p opula tion of zero-age massive binaries. 
Similarly, the conclusions of iPinsonneault fc Stanekl (120061 ) are based upon a small number 
of eclipsing systems in the SMC, and are, therefore, most sensitiv e to systems of equal mass 
and radius so that no firm conclusions can be drawn (iLucyl 120061 ) . 



In this pape r, we analyze data fro m the radial velocity survey of 120 O and early-B stars 
in Cygnus OB2 (IKiminki et al.l 120071 ) to help resolve these disparate conclusions regarding 
the bi nary characterist i cs of massive stars. The new data more than double the sample size 
of the IGarmany et al.l (Il980l ) study, include massive stars down to early B spectral types, 
and often involve more epochs of observations. Although orbital parameters for a large 
number of systems are not yet available, we use the observed distribution of radial velocities, 
in conjunction with Monte C arlo models, to test seyeral popular formulat i ons of F, f(q), 



and f (r) from the literature (IPinsonneault fc Stanekl 120061 ; lHogeveenlll992l ; IGarmany et al. 
19801 ). Using only the raw velocities of a large sample of stars from one young OB association 
mitigates many of the selection and evolutiona ry effects r esponsible for conflicting results 
from more general binary surveys and catalogs (IToutl 1 1 99 ll ). 



We find that we are able to constrain the probable values for a and (3 to a narrow region 
of parameter space, subject to reasonable assumptions about the boundary conditions of the 



mass ratios and orbital separations. Only select combinations of a, /3, 



^out 



and F can 



reproduce the levels of velocity variation seen in the data. Some canonical parameterizations 
for binary characteristics are inconsistent with the data and can be ruled out. Our immediate 
intent is to provide improved input parameters for population synthesis models that predict 
the rates of energetic phenomena in the Milky Way and at cosmic distances. We also expect 
that these results will have implications for the formation scenarios of massive stars which 
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are c urrently under debate (jBonnell et al.l Il998l ; IWolfire fc Cassinellil 119871 ; iKrumholz et al 
2005h . 



The Data 



Kiminki et al.l (120071 ) report on a 6-year radial velocity survey of 120 early-type probable 
members of the Cygnus OB2 Association. Cygnus OB2 was chosen as a target for long term 
radial velocity study to determine the binary fraction and mass ratio distribution among 
massive binaries in a populous young OB association where the stars are of comparable age 
and stellar evolutionary effects are minimized. However, even in the ~2-3 Myr old Cygnus 
OB2, the pr esence of a variety of evolved stars sugges ts that the star formation process was 
non-coeval (IMassey fc Thompson! Il990l ; lHansonll2003l ). 



Spectroscopic observations from the Lick Observatory Shane telescope, the Keck Ob- 
servatory Keck I telescop^, the Kitt Peak National Observatory^ WIYN0 telescope, and 
Wyoming Infrared Observatory 2.3 m telescope (WIRO) over the period 1999 July to 2005 
October were used to measure radial velocity variations of stars earlier than spectral type 
B3. Sample stars were observed on at least 3 epochs over this period. Some stars have 
as many as 18 epochs of data. The mean and median number of epochs are 7.6 and 7.0 
respectively. The time sampling is irregular, as imposed by telescope schedules and weather 
conditions. Sampling is additionally modulated by the observing season for Cygnus (June 
through Novem ber). The data were r educed and analyzed using cross-correlation techniques, 
as described by iKiminki et al.l (120071 ) . 



Figure [T] shows the analytical two-body relationship between velocity semi-amplitude 
and orbital period for binary systems with a 10 M Q primary. Solid lines indicate the loci 
of systems with equivalent mass ratios, q = M2/Mi=1.0, 0.1, 0.01. Dotted lines indicate 
loci of common orbital separations, r = 0.1, 1.0, 10 AU . The dark gray area in the upper- 



1 Some of the data presented herein were obtained at the W.M. Keck Observatory, which is operated 
as a scientific partnership among the California Institute of Technology, the University of California and 
the National Aeronautics and Space Administration. The Observatory was made possible by the generous 
financial support of the W.M. Keck Foundation. 

2 NO AO is the national center for ground-based nighttime astronomy in the United States and is operated 
by the Association of Universities for Research in Astronomy (AURA), Inc., under cooperative agreement 
with the National Science Foundation. 

3 The WIYN Observatory is a joint facility of the University of Wisconsin, Madison, Indiana University, 
Yale University, and the National Optical Astronomy Observatories. 
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right demarcates the prohibited region of parameter space where q > 1. The light gray area 
indicates the portion of parameter space sampled in this survey. The data are sensitive to 
primary velocity semi-amplitudes > 5 — 10 km s _1 , periods shorter than ~ 2 x 6 years and 
orbital separations <~ 2 x 10 AU. 



The data presented in lKiminki et al.l (120071 ) show at least several dozen stars with strong 
radial velocity variations at the level of > 20 km s -1 and dozens of additional velocity- 
variable candidates with semi-amplitudes down to the survey sensitivity of ~ 5 — 10 km s _1 . 
There are also ~ 6 probable double-lined spectroscopic binaries, but the individual compo- 
nents are poorly resolved in most of our data, so these systems are, for present purposes, 
treated as single-lin ed binaries. Fo r this analysis, we removed from the main sample of 120 
objects (Table 5 of iKiminki et al.l (120071 )) six stars later than B4. The remaining sample 
used for analysis is 114 obj ects. These include the stars identified as probable members by 
Massev fc Thompson! (|1990r) pl us additional earl y B type memb e rs ide ntified photometrically 
by iMassey fc Thompson! (Il990l ) and classified by IKiminki et al.l (120071 ) . Given the significant 
uncertainties in the photometric distances, we have undoubtedly included a few foreground 



and b ackground stars. The distribution of photometric distances (Figure 5 of IKiminki et al. 



(120071 )) is nearly Gaussian with only a couple of outliers, suggesting that the sample does not 
suffer significant contamination from field stars. Including a few field objects in the 114-star 
sample will not measurably impact the conclusions, even if the binary characteristics of field 
OB stars are different than the Cygnus OB2 members. 

Figure [2]plots the log of the assigned spectroscopic masses (IMartins et al. hood : iHumphreys Sz McElroyi 

1984 ) for each primary star versus the quantity Vh = 0.5|V^, aa; — V m i n \, where V max and V m i n 



are the maximum and minimum observed velocity for a given star. Vh is a measure of the 
velocity semi-amplitude of the primary, albeit an imperfect one because the velocity curves 
are not generally well sampled at all phases. We regard Vh as a lower limit on the true 
projected velocity semi- amplitude. Filled circles in Figure [2] denote main sequence OB stars 
while triangles denote evolved stars. The error bars show the mean uncertainty in the veloc- 
ity measurements over all epochs for a given star. A small dispersion has been added to the 
masses in a few hea yily populated mass bins to improve clarity. One 03If star (MT457 in 
the nomenclature of Massey fc Thompson ( 199ol )) at M = 80 M falls outside the maximum 
plot range. The typical velocity precision is ~ 5 — 15 km s -1 , and it varies somewhat with 
source luminosity and observatory/instrument combination used. For example, the radial 
velocity uncertainties obtained with the WIYN telescope + Hydra spectrograph are often 
better than those obtained with the Lick 3 m telescope + Hamilton echelle spectr o graph 
because of the higher signal to noise ratio achievable with the former. IKiminki et al.l (120071 ) 
provide a more complete description of these data and place a lower limit on the binary 
fraction at 30% - 42%. 
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Figure [3] shows a histogram of the observed velocities and uncertainties. This Figure 
illustrates the distribution of velocity dispersions, V rms , calculated from the multiple mea- 
surements of the 114 OB stars (solid line) along with the distribution of mean velocity 
uncertainties, a v (dashed line). The dotted line shows the distribution of Vh- The lowest 
velocity bin from to 5 km s _1 is sparsely populated because observational errors scatter the 
data into higher velocity bins. The maximum observed semi-amplitudes fall mostly between 
10 and 40 km s _1 , with a significant tail toward higher velocities out to ~90 km s _1 . The 
uncertainties lie in the characteristic range 5-15 km s _1 . 

For comparison, Figure H] shows the distribution of properties for O-type systems from 
Garmany et al.l (Il980l ). The left panel shows a histogram of the velo city semi-amplitudes fo r 
the new observations of O stars presented in Table 3 and Figure 1 of iGarmany et al.l (ll980l ). 
Note the differences compared to our Figure [3J Most striking in the Garmany et al. sample 
is the relative lack of stars in the velocity range 20-80 km s _1 and the relative excess of 
star s in the range 8 -120 km s _1 . Significant selection affects may drive these differences. 
The IGarmany et al.l (jl980j) sample was magnitude limited and contained stars at a variety 
of evolutionary stages, including a large fraction of supergiants. For such systems, mass 
exchange and transfer of angular momentum between orbital and rotational c omponents may 



have a ltered the original binary mass ratios and separations. Therefore, the IGarmany et al. 



(119801 ) sample may not be representative of the general population of massive binaries at 
zero age. The middle pa nel of Figure HI i s a hi stogram of mass ratios, q, for the O type 
binaries from Table 3 of IGarmany et al.l Jl980h fl The histogram shows a preference for 
mass ratios near unity and is consistent with a flat (a = 0) or rising distribution (a > 0) 
with increasing mass ratio. The right panel of Figure H] shows a derived semi-m ajor axis 
distribution for the O type binary systems from Table 3 of IGarmany et al.l (Il980l ) assuming 
an average orbital inclination of 60°. The dashed line is a best fit Gaussian curve which 
peaks near log(r) = —0.8 or ~ 0.15 AU. This panel shows that the O-type binaries from 
Garmany et al. are preferentially th ose with small orbital s eparations, as might be expected 
since these are the easiest to detect. IGarmany et al.l (Il980l ) did consider whether the lack of 
systems with small velocity amplitudes implied a true absence of such systems or a merely 
limitation on their detection. They concluded that such low-amplitude systems should have 
been detected as part of their survey and that the preference for equal mass components was 
a real effect. 



4 Hcre, we have included only the stars with new observations denoted by an asterisk in their Table 3. 
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3. Modeling Binary Characteristics 
3.1. Monte Carlo Modeling of Radial Velocities 



We wrote a Monte Carlo code to simulate the radial velocity variations expected for 
Cyg OB2 primaries based upon analytic prescriptions of binary mass ratios and orbital 
separations. We initially assumed a binary fraction of F = 1.0 so that every primary star 
has one companion. Statistically, this means that single stars are treated as systems with 
very low mass companions, very large separations, or very low inclinations. However, this 
is a u seful approach to ob tain an initial lower limit on a. For each of the 114 primary 
stars (IKiminki et ajj 120071 ) of mass Mi and observed velocity variations characterized by 
Vh and V rms , the Monte-Carlo code randomly assign a secondary star of mass M 2 based 
upon a distribution of mass ratios, f(q) where q = M 2 /Mi. The mass ratios are drawn 
from a population described by f(q) oc q a where a is the power-law index of the mass ratio 
distribution over the range qi ow < q < 1. Initially, we adopt qi ow = 0.02 so that the lowest 
mass secondaries for our least-massive ~8 M primaries are 0.16 M Q and ~1 M & for the 
handful of very massive 50 M & stars. Such masses are sufficiently low that the secondaries 
produce negligible radial reflex motion in their primaries — well below the velocity sensitivity 
of the data. We consider values of a between -3.0 and 2.0. a > describes a secondary star 
mass distribution peaked toward the primary mass, a = corresponds to a flat distribution 
of secondary masses. The Monte-Carlo code assigns an orbital separation, r, drawn from a 
distribution described by /[log(r)] oc [log(r/rj n )]^ with rj n < r < r out and r in A.U. Opik's 
Law (/? = 0) corresponds to a parent population of orbital separations uniformly distributed 
in log(r) between log(r in ) and log{r out ). [3 = ±3 corresponds to distributions peaked toward 
large/small separations, respectively. 

Because the power-law slope, (3, that best describes the orbital separations of the data 
will depend upon the chosen boundaries, r in and r out , we explore a small 3x3 grid of values: 
r in =0.020, 0.063, & 0.200 AU and r out =100, 1000, & 10000 AU. These inner limits 
correspond to 4.2, 13, and 43 R Q . For a main-sequence B0 star, representative of our targets, 
these inner limits correspond to 0.56, 1.7, and 5.8 R*. Thereby, the smallest two values for ri n 
are firm lower limits, representing separations where significant mass transfer and common 
envelope evolution would take place. Although these values are physically implausible in 
general, we include them to illustrate the sensitivity of the simulations to this parameter. 
The third value, 0.20 AU, is the largest reasonable choice for rj n , as it i s larger than several o f 
the orbital separations for short-period Cyg OB2 binaries reported in IKiminki et al.l (120071 ) 
and lGarmany et al.l (119801 ) (~ 0.15 AU). The minimum adopted outer limit, r out = 100 AU is 
small compared to o bserved orbital separation s in wide binaries which may extend to several 



thousand AU (e.g., iKouwenhoven et al.l 120051 ; iDuquennoy fc Mayor! Il99ll ). The maximum 
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adopted outer limit, r out = 10000 AU, is sufficiently large that few systems are likely to 
remain bo und at such distances , espec ially in the dense stellar environments where OB stars 



are born. iLepine fc Bongiornol (120071 ) report that among nearby Hipparchos binaries, the 
Opik's Law distribution of separations fits the data out to a maximum of r = 4000 AU, 
beyond which the number of companions drops, consistent with gravitational disruption. 

In the interest of minimizing free parameters, we assume that the orbits have zero ec- 
centricity. The resulting period and velocity semi-amplitude of the simulated binary system 
is computed N iter times, where N iter is typically 100 Monte-Carlo iterations. For each sim- 
ulated system, the velocity curve is sampled at a random phase angle N f, s times, where 
3 < N b s < 18, the actual number of observations of a given star. The code assigns a random 
inclination angle, i, for each system, where i is generated by allowing the angular momentum 
vector of the simulated system to lie anywhere on a sphere. The mean observational velocity 
uncertainty for each individual star, normally distributed about zero, is then added to the 
simulated velocity. As a test, we also performed simulations using the actual uncertainties 
for each measurement of each star and we found that the results differed insignificantly from 
using the average uncertainty. 

For each Monte Carlo run with a given pair of a and (3, the simulated distribution 
of Vh and V rms for 114 stars is compared to the observed distribution using a two-sided 
K-S test. The probabilities that the simulated and observed distributions of Vh and V rms 
(P(Vh) and P(V rrns ), respectively) are drawn from the same population are averaged over 
100 Monte Carlo iterations and tabulated for later analysis. We found that both velocity 
metrics, Vh and V rms , yield similar results and provide similar probabilistic constraints on 
the binary characteristics. While Vh is a better measure of the true semi-amplitude, it is also 
less robust than V rms against outliers or measurements with large uncertainties. Therefore, 
in subsequent discussion and figures we use the average probability, P a , obtained from the 
arithmetic mean of P{Vh) and P(V rms ). 

In Section 4 we will discuss the contraints imposed by this modeling appraoch. The 
Figures and discussion will show that there is degeneracy between a and (3 — in general, 
large a coupled with large (3 reproduce the velocity data approximately as well as small a 
coupled with small (3. Additional constraints are needed. One such constraint we consider is 
that Opik's Law ((3 = 0) obtains. Another is the fraction of type Ib/c supernovae predicted 
by each model. 
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3.2. Monte Carlo Modeling of Type Ib/c Supernva Rates 



The massive binary population can also be constrained by observations of energetic phe- 
nomenae that arise from massive binary progenitors. A well-studied example is the fraction 
of core-collapse su pernovae that o ccur i n stars that have lost their hydrogen envelopes — type 
Ib/c supernovae. iPodsiadlowskil (119921 ) argued that most Ib/c supernovae are produced in 
binary star systems where mass transfe r removes most /all of the hydrogen envelope, unveil- 
ing a type Ib/c supernova progenitor. Podsiadlowski ( 19921 ) assumed that 37% of binaries 
were close binaries that would undergo some mass transfer. 

The fraction of type Ib/c supernovae predicted by population sysnthesis models must 
agree with observations and must self-consistently produce the levels of velocity variation seen 
in massive primaries, such as those in Cyg OB2. The fraction of type Ib/c supernovae rel ative 



to all core-collapse sup ernovae is fn,/ c = N Ib / c /(N Ib / c + Nu). ICappellaro et al.l (119991 ) and 
Mannucci et al.l (120051 ) give the observed fraction of type Ib/c supern ovae in Sbc/d spiral 



galaxies such as the Milky Way as 15%. The more recent works of iLi et al.l (120071 ) and 
Leaman et al.l (120071 ) revise this estimate upward to 30%±11%. Using binary population 
synthesis calculations, we can apply these fractions to further constrain the parameter space 
in our study. Any part of parameter space that produces a fraction of type Ib/c supernovae 
inconsistent with the observed fraction can be ruled out. 



We used the binary synthesis code developed by iFryer et al.l (119981 . Il999l ) to calculate 
the total core-collapse supernova rate and the number of type Ib/c (stars that have lost 
their hydrogen envelopes at collapse) supernovae. Like Podsiadlowski (1992), we find that 
the primary formation path for type Ib/c supernovae begins with a close binary where the 
more massive primary envelops its companion in a common envelope phase. The companion 
spirals into the envelope of the massive primary, ejecting the hydrogen envelope. After the 
common envelope phase, the primary will continue to lose mass in a wind. For this paper, 
we will assume any star that loses its entire hydrogen envelope has the potential to become a 
type Ib/c supernova. Because mass is lost in winds as well, even systems (especially those at 
high metallicity) that merge in the common envelope phase without ejecting the hydrogen 
envelope can make type Ib/c supernovae. Systems with more extreme mass ratios (i.e., 
q < 0.2) are more likely to go through a common envelope phase and systems with such 
low-mass companions dominate the formation rate of type Ib/c supernovae. Figure O shows 
the range of primary masses and orbital separations that form Ib/c supernovae for 4 separate 
values of the system mass ratio. As q increases, the system tends to undergo more stable 
mass transfer without a common envelope phase. Most of these systems do not become type 
Ib/c supernovae. These simulations are based on our standar d wind assumptions in which 
wind mass loss rates are set to 10% of the lWoosley et al.l (119951 ) values. Beyond 32 M Q , these 
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primary star winds will also produce type Ib/c progenitors, regardless of companion mass. 
However, this corresponds to only a small fraction of our type Ib/c supernova rate. 

The other requirement for a star to produce a type Ib/c supernova is that it actually 
explode as a supernova. We assume any primary with an initial mass above 9 M will 
produce a supernova. Although the total rate of type Ib/c supernova does depend on this 
choice, the dependence on the type II rate is nearly identical and varying this value in the 
8 — 10 Mq range leads to variations in the rate of 2-3% for the bulk of our parameter study. 
We also assume that stars that have helium cores in excess of 11M Q collapse directly to 
black holes and don't form supernovae. Varying this value of the helium core mass over the 
8 — 15 M range yields variations in the rate of less than 1%. This is not surprising because 
the number of such high-mass systems is quite small. 

In this paper, we focus our binary population synthesis studies on the mass distributions 
of the stellar components and the orbital separation. These distributions are varied accord- 
ing to range of parameters studied in the analysis of our observations and constitute over 
10,000 different population synthesis calculations each including 100,000 binaries for each 
calculation. In all of these calculations, we assume the eccentricity is zero. If we instead 
assume the eccentricity follows a flat distribution between and 0.8, we can increase the 
rate by 5%. 

We have already discussed the dependence of our ratios on some of the many free pa- 
rameters in population synthesis: e.g. critical mass for supernova formation. Let's briefly 
discuss some of the remaining free parameters. We assume a single Maxwellian kick distri- 
bution with a mean value of 300 km s -1 . Because nearly all of the type Ib/c supernovae are 
formed in hydrogen envelope stripping in the common envelope inspiral of the companion 
star into the massive primary, the rate does not depend strongly on the choice of this kick. 
We have varied the value of this mean value from 100 km s -1 to 600 km s -1 as well as using 
the bimodal kick distribution of Fryer et al. (1998) and the rate does not change by more 
than 5%. For mass transfer, we use the same formalism as Fryer et al. (1999) with mass 
transfer parameters q;mt = 1-0 and Pmt = 0.5 and a common envelope efficiency «ce = 0.5. 
Varying these values for a wide range of values (e.g. 0.2 < a CE < 1.0) yields variations in 
the ratio of 10-30%. 

The final two uncertainties studied in our calculations are our assumed stellar radii and 
stellar mass loss rates. Here again, we use the standard Fryer et al. (1998) values. But if 
we allow a factor of 2 variation (both increase and decrease) in these values, the ratio of 
type Ib/c supernovae to type II supernovae can vary by nearly 50uncertainty is one of the 
dominant uncertainties in the population synthesis of binary systems: the error in the stellar 
radius. If the stellar radii are smaller, we produce fewer common envelope systems with a 
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larger fraction of the remaining common envelope systems leading to mergers prior to the 
ejection of the hydrogen envelope. This error dominates our 50% error. Clearly we can vary 
the population synthesis parameters to make many of the models fit to the data. But if we 
assume we know the stellar radii and take standard values for common envelope efficiencies, 
our calculation errors can be limited to 5-10% effects, allowing us to place strong constraints 
on the assumed distributions of initial masses and separations. 



4. Results of Modeling 

4.1. Binary Characteristics for single power-law distributions of q and r 

Panels within Figure M show probability contours (solid lines) depicting the likelihood, 
P a , that a given combination of a and (3 reproduces the distribution of velocities in the 
data for a minimum mass ratio of q — 0.02 and binary fraction of F — 1.0. The 3x3 
grid corresponds to r in = 0.020 AU (left column), r in = 0.063 AU (middle column), r in = 
0.200 AU (right column), and r out = 10000 AU (upper row), r out = 1000 AU (middle row), 
r out = 100 AU (lower row). The horizontal line in each panel marks the nominal (3 = 
(Opik's Law) distribution of separations. The dotted contours in each panel depict the 
predicted fraction of type Ib/c supernovae in our Monte Carlo models. These contours run 
from 10% to 40% in increments of 5%. 

The panels in Figure M demonstrate the sensitivity of the best-fitting solutions for a 
and (3 to the adopted r in and r out . Each panel contains a crescent-shaped ridge of highest 
likelihood contours (P a > 50%). Values outside of this region drop to probabilities < 10% 
over much of the plotted parameter space. These contours show that the best-fitting a 
and (3 are correlated — larger orbital separations are required (i.e., larger (3) when mass 
ratios are peaked toward unity, while smaller orbital separations are required for mass ratio 
distributions with secondaries drawn from a field star IMF. The ridge of peak likelihood 
separates parameter space into two regimes: to the left of the ridge of solid contours the 
models fail to produce enough high-velocity systems and are inconsistent with the data; to 
the right of the peak contour ridge, the models produce too many high-velocity systems 
relative to the data. 

For a flat distribution of orbital separations described by Opik's law, the best-fitting 
values for a lie in the broad range —1.2 < a < 1.0 depending on the choice of radial limits. 
For the most reasonable choices r in =0.063-0.200 AU and r out = 1 000 AU, a lies in the range 



0.6 - 0.0, broadly consistent with obser vational estimates of (IKouwenhoven et al.l 12005 



Kouwenhoven 1120061 ; iGarmany et al.lll980l ). The dashed contours depicting fn,/ c =0.25-0.35 
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intersect the (3 = line and the solid co ntour ridge, indicating a Ib/c fraction about t wice 



the older observational estimates of 15% (IMannucci et al.l 12005 



Cappe 



aro et al. 



1999) but 



consi stent with the most recent determinations of 30%±11% (iLi et al.l 120071 ; iLeaman et al. 
20071 ). Ib/c fractions of fjy c <0.20 are inconsistent with the data in all panels except the 



upper-right which shows the (perhaps unrealistically) large radial limits of ri n = 0.20 AU 
and r out = 10, 000 AU. It is noteworthy that, in all of the panels, mass ratio distributions 
described by a < —2 are unlikely (contours levels <30%) even if the distribution of orbital 
separations favors close companions (small (3). Furthermore, any < a < — 1 for the most 
plausible radial limits in the middle or upper middle panels would require /n,/ c > 0.40 and 
(3 < —0.3, in disagreement with observational results. 

Figure [7] shows the distribution of observed and simulated velocity semi-amplitudes, Vh, 
(solid and dashed lines, respectively) for a subset of models that produce P a ~ 60% from 
the contour ridge in the middle- right panel of Figure El The left panel of Figure [7] shows a 
simulation with a = 0.4, (3 = 0.10, and the right panel shows a = —0.6, (3 = —0.20. Each 
panel is labeled with the probability that the two histograms are drawn from the same parent 
population. Both combinations of a, (3 yield probabilities of P(Vh) — 0.60 and lie along the 
ridge line of peak likelihood. There is good agreement between the data and simulations 
based on the two-sided K-S "D" statistic which quantifies the maximum difference between 
the cumulative distribution functions of the two data sets. An analysis of the cumulative 
distribution functions shows that, for these best-fitting models, the maximum differences 
occur in the range 10-40 km s -1 and 80-100 km s _1 in roughly equal proportion. This 
means that, in some Monte Carlo iterations, differences in the high-amplitude tail of the 
velocity distributions limit the agreement between the models and data; in other iterations, 
the models fail to produce the observed velocity distribution on the low-amplitude side 
of the histogram. The lack of systematic trends in the velocity at which the K-S "D" 
maximum occurs may be taken as an indication that the comparison is limited by sample 
size and measurement uncertainties. If, for example, the K-S "D" maximum difference always 
occurred at small Vh amplitudes, we would suspect a systematic problem with the adopted 
velocity uncertainties which are responsible for the objects in the lowest velocity bins. If, on 
the other hand, the K-S "D" maximum difference always occurred at large Vh amplitudes, 
we might suspect the presence of a second population of close binaries which was responsible 
for the large-amplitude systems. We will return to this latter possibility in a subsequent 
section. 

Figure [S] shows contour plot comparisons of the data with Monte Carlo simulations 
in the same manner as Figure EJ except for a binary fraction of F = 0.8. The overall 
level of agreement between the models and the data is lower in each panel by ~10% when 
compared with the corresponding panel of Figure [6j The contours are qualitatively similar 
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to those in Figure [6] in that the ridge of peak likelihood reaches P a ~ 50% in most panels, 
but it is shifted slightly toward larger a and smaller (3. Variations between panels once 
again highlight the sensitivity of the best-fitting solutions to the adopted radial boundary 
conditions. For instance, the middle-right panel (r in = 0.20, r out = 1000) shows a peak 
likelihood of P a > 60% for an Opik's Law distribution of r, but only for a > 0.5, meaning 
that mass ratios must be peaked toward unity. This is because, with fewer binaries, the 
observed levels of velocity vari ation can only be repr oduced with larger mass ratios and/or 
smaller separations. Thus, the iGarmany et al.l (119801 ) mass ratio distribution peaked toward 
unity can be consistent with the data, but only if the binary fraction is F < 0.8 and only for 
somewhat large value s of = 0.20 AU, i.e., larger than many of the systems observed by 
Garmany et al.l (119800 . For the most plausible radial limits (rj n = 0.063, r out = 1000) in the 
center panel, the most probable value for a is ~ 0.0 for (3 = 0. This implies frh/r. =0 .20-0.25, 



at the lower end of current observational constraints (ILi et al.l 120071 ; iLeaman et all 120071 ). 
Any of the other panels for plausible radial limits (upper-middle or middle-right panels) 
require fjb/ c <0.20, a value that is increasingly inconsistent with observational limits. 

Figure [9] shows contour plot comparisons of the data with Monte Carlo simulations in the 
same manner as Figure for binary fractions of F = 0.6. Compared with the corresponding 
panels of Figure [8], the solid contours in each panel shows lower levels of agreement between 
models and data by about 10%. This is especially true in the middle and upper-middle 
panels which represent the most plausible radial boundary limits (rj n = 0.063 — 0.200, 
r out — 1000 — 10000). These panels also show that mass ratio distributions with a < — 2 
are highly inconsistent with the data even if the the orbital separations are preferentially 
very small {(3 < —1) and rj n = 0.020 AU. Thereby, the field IMF distribution for secondary 
star masses is highly improbable given that Ti n = 0.020 would require many systems in 
common envelopes even for main-sequence primaries. Reasonable agreement with the data 
(P a > 50%) is only achievable with a > and with r out = 100 AU (lower row). Because 
r out — 100 AU is implausible (no long-period companions!) and any of the panels for 
the most probable radial limits (middle or upper-mid dle) require fib/ c <0.15, at odds with 
current estimates (ILi et al.l 120071 ; ILeaman et all 120071 ). we conclude that Figure [9] contains 
no acceptable portion of parameter space where a good fit between the models and data can 
be achieved. A binary fraction as of F = 0.6 or lower is highly improbable in Cyg OB2. 

In summary, the sequence of Figures El [HI El shows the degeneracies between a and 
(3 and the sensitivity of the best-fitting solutions to the adopted boundaries of the power- 
law parameterization of separations, f(r). Distributions of q favoring extreme mass ratios 
a < —2 universally require very small separations (J3 < —1) in disagreement with observa- 
tions supporting the canonical flat Opik's Law distribution. For plausible inner and outer 
boundaries, the most probable values for a lies in the range —0.6 - 0.0 for F = 1.0 and 0.0 - 
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1.0 for F = 0.8. Binary fraction F < 0.8 would imply low fractions of type Ib/c supernovae, 
fib/c <0.20, inconsistent with current estimates. 



4.2. Comparison with the Hogeveen Distribution of q 



Hogeveenl (11990l . Il992l ) proposed that /(g) oc g~ 2,0 2 7 for q > g where g = 0.3 for B 
star primaries, and f(q) oc q° (i.e., flat) for q < q . The relative numbers of systems above 
and below go is not given. We used an adap ted Monte Carlo code to perform comparisons 
between the data and the predictions of the iHogeveenl (119921 ) mass ratio distribution. We 
adopted a two-component power-law distribution of q with variable normalization for each 
component so that they yield the same /(g) at q . The component with q < q has fixed 
a(q < go) = while a(q > go) is allowed to vary freely. This effectively means that the 
fraction of systems with g > g varies with a. For a = 2.0, 90% of systems have g > go; 
for a = 0.0, 70% of systems have g > g ; for a = —2.0, 48% of systems have g > g . We 
compared the velocity distribution of the data to the expectations of Monte Carlo simulations 
in the same manner as for the simple power-law parameterizations of f(q ) and f(r) above . 
We note that in the limiting case a(q > go) = 0, the predictions of the IHogeveenl (119921 ) 
mass ratio distribution and the power-law formulation discussed above are identical, within 
the noise of the Monte Carlo procedure employed. 

Figure [TUl shows probability contour plots as a function of a(q > g ) and (3 for a range 
of ri n and r^t as in Figure [61 The ridge of peak likelihood exceeds 60% in some panels. 
In particular, the center (rj n = 0.063) and middle-right (rj„ = 0.20) panels show best 
agreement with the data. Both panels have r out = 1000 AU . The center panel shows 
that for the nominal Opik's Law distribution of separations, the most probable values for 
a lie in the range —3.0 < a < —2.0, consistent with the Hogeveen formulation for q > go. 
However, the predicte d fraction of type Ib/c supernovae is //&/ c =0.35-0.40, at the upper edge 
of current estimates (ILi et al.ll2007l ; iLeaman et al.ll2007l ). In most other panels, a < —2 is 
inconsistent with the data. The middle-right panel shows that for the Opik's Law distribution 
of separations, the most probable values for a lie in the range —1.0 < a < 0.0. Hence, the 
most probable value for a for the Hogeveen mass ratio distribution is much more sensitive 
to the choice of r« n than for the simple power-law distribution of q. This is because the 
peak contour ridge in each panel is nearly horizontal, a consequence of the fact that the 
actual distribution of mass ratios over the whole range q =0.02 - 1 varies comparatively 
little because the portion with q < go is unchanging. Under the Hogeveen formulation, the 
two-component power law with a(q < go) = and a(q > go) = —2 essentially mimics a single 
power-law of intermediate slope, a ~ —0.6. The predicted fractions of type Ib/c supernovae 
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in the middle and middle-right panels are fib/ c =0.35-0.40 and fjj,/ c ~0.30 respectively. 

Figure [TT] shows probability contour plots as a function of a and (3 for a range of r in and 
r out as in Figure [TD], except for a binary fraction of F — 0.8. The overall level of agreement 
between the models and the data is lower in each panel by ~ 10% when compared with the 
corresponding panel of Figure [10J The contours of peak likelihood are shifted toward larger 
a and smaller (3 relative to Figure [TOj In the center panel, the peak probabilities exceed 
50%, and for an Opik's Law distribution o f separations, the best fitting value of a lies in the 
range —1 < a < 0.3, inconsistent with the iHogeveenl (119921 ) value of a = —2. In the middle- 
right panel, the peak probabilities exceed 60%, but Opik's Law would require a > 0.8. The 
predicted fraction of type Ib/c supernovae is fjb/ c =0.20-0.25, at the lower end of current 
estimates. The lower-middle panel shows probabilities exceeding 60% for Opik's Law and 
a ~ —2. However, this requires r out = 100 AU, and we argue previously that such a small 
outer bound on companion separations is not supported by observations. In summary, the 



Hogeveenl (119921 ) mass ratios distribution and a binary fraction of F — 0.8 produces slightly 



poorer overall agreement with the data. The best fitting models consistent with Opik's Law 
require — 1 < a < 1 and r in = 0.063-0.200 and r out ~ 1000 AU. Simulations with F < 0.6 
show comparatively poor agreement with the data, similar to Figure [HI and would require 
fib/c <0.15, a value inconsistent with current measurements. 



4.3. Comparison with a "Twin" Component Model Distribution of q 



Given ample evidence in the literature for a significan t "twi n " population among short 
perio d binaries (jPinsonneault fc Stanekl 120061 ; iTokovininl |2000| ; iLucy fc Riccd Il979l ; iLucy 
20061 ). we used a modified Monte Carlo code to simulate the expected velocity amplitudes 
from a two-component q distribution. The dominant component has a power-law distri- 
bution, f(q) oc q a and /(logr) oc [log (r/r in ) Y between the limits r in = 0.063 AU and 
r out = 1000 AU, with a and (3 being free parameters, as described previously. The twin 
component comprises a fraction, Tf rac , of the systems. We explore a range of twin fractions, 
Tfrac =0.40, 0.25, and 0.10. Allowing for the possibility that F may be less than unity, we 
define Tf rac to be the fraction of systems with q ~ 1 among those that are binary systems. 
Motivated by observations that "twin" systems have small separations, presumably because 
they form within a causally connected region, we adopt a log-normal distribution with mean 
separation, < r Twin >, and Gaussian width, o"i ogr . We explore a range of typical separations, 
< rTwin >=1.0, 3.9, and 15 AU (i.e., < logrxwin >=0.0, 0.6, and 1.2 ). The smallest of these 
mean separations, < rTwin >= 1 AU, is representative of the orbital distances required to 
produce OB star twins with orbital periods in the range few-60 days like the samples used 
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by (jPinsonneault fc Stanekll2006l ; iTokovininl I200CH : iLucy &: Riccdll979l ) to infer the presence 



of a twin population. We adopt (J\ ogr = 0.6 for the width of the twin population's orbital 
separation distribution. (Exploration of larger and smaller values for a\ ogr yielded poorer 
levels of agreement with the data and we do not show these results in the interest of brevity.) 
The mass ratios for the twin component are uniformly distributed between q =0.95-1.0. The 
lower mass limit for the power-law component described by f(q) oc q a is qi ow = 0.02. 

Figure d2] shows a 3x3 grid of contour probability plots T) rac =0.40 (left column), 7} rac =0.25 
(middle column), Tj mc =0.10 (right column), and < r Twin >=15 AU (upper row), < r Twin >=3.9 
AU (middle row), < rTwin >=1.0 AU (lower row). The horizontal line in each panel de- 
notes Opik's Law (fl = 0). The panels in the right column having Tf rac = 0.10 closely 
resembles the center panel of Figure [6] where no twin component is included, indicating that 
twin fractions of <10% do not measurably alter the results. In general, the lower row with 
< fTwin >=1.0 AU shows relatively poor agreement with the data — small mean separations 
yield too many large-amplitude systems. The peak probabilities exceed 60% in several pan- 
els, most significantly in the middle panel corresponding to < rTwin >=3.9 AU, T/ rac =0.25 
and in the middle-left panel corresponding to < rTwin >=3.9 AU, Ty rac =0.40. In these 
panels Opik's Law specifies best-fitting values for a in the range —3.0 < a < —1.3. The 
predicted fractions of type Ib/c supernovae are large in both cases, exceeding 0.40. In the 
middle-left panel showing < rTwin >=3.9 AU, Tf rac = 0.40, there is a broad region where 
the peak probabilities exceed 60% and reach 70%. In this case, if Opik's Law obtains, then 
a < —2, consistent with the expectations of randomly pairing secondaries from the field 
IMF with massive primaries. Figure d2] illustrates that two-component q prescription of 



Pinsonneault fc Stanekl (120061 ) (45% of systems with q > 0.95, the remainder with a = 0) 
is not consistent with the data unless (3 > 0.5. That is, the probability of a = 0, /3 = 
is < 30% for all of the panels. The presence of any significant twin component appears to 
require a < — 1. 

Figure [13] shows a 3x3 grid of contour probability plots similar to Figure [12] but for 
a binary fraction of F — 0.8. The peak contours once again exceed 60%, primarily in 
the upper two rows of panels. The ridge of peak likelihood is shirted toward larger a and 
smaller (3 relative to Figure [T2] In most panels, the region with P a > 0.60 is smaller than 
in the corresponding panels of Figure [12] indicating slightly poorer overall agreement with 
the data. For an Opik's Law distribution of separations, the most probable values for a 
lie in the range -1.5 - 0.0 for twin fractions Tf rac =0.25-0.40 and < rTwin >=3.9— 15 AU. 
Models with < r Tw i n >=1 AU produce generally poorer ag reement with the data, having 



peak likelihood < 50%. The two-component q prescription of iPinsonneault fc Stanekl ((2006) 
(45% of systems with q > 0.95, the remainder with a = 0) is not consistent with the data 
unless < rTwin >— 15 AU. At such large distances, OB star twins would have much longer 
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periods (~ 10 yr) than the ~30-day systems us e d to infer the presenc e of a twin population 
( IPinsonneault &: Stanek! 1200a : iTokovininl l2000l : iGarmany et al.l Il980l ) . For (3 = and in 



panels with P a > 0.60 (middle, middle-left and upper-left), the predicted fr actions of type 
Ib/c supernoyae are fjb/ c =0.20-0.30, in broad agreement with observations (ILi et al.l 120071 : 
Leaman et al.ll2007l ). 



Figure [14] shows a 3x3 grid of contour probability plots similar to Figure [12], but for a 
binary fraction of F = 0.6. The ridge of peak likelihood is now shifted ever farther toward 
large a, and the peak probabilities do not much exceed 50% in most panels, indicating poorer 
agreement with the data compared to Figure [131 In those panels (mi ddle, middle-left, upper- 
middl e) , Opik's Law would require a > 0, in disagreement with the IPinsonneault fc Stanek 
(120061 ) and field-star IMF distribution. Furthermore, the type Ib/c supernova fraction is 
fib/c <0.20, in disagreement with observations. Binary fractions F < 0.6 appear improbable. 



4.4. Comparison with a Miller-Scalo Secondary Mass Distribution 



Power-law parameterizations of the stellar initial mass function are appropriate only over 
a selected mass range. More realistic descriptions of the IMF include multiple power-law 



components with slopes t hat flatten toward lower masses (e.g., iMiller fc Scald Il979l ; IScalo 
198a : iKroupa et al.lll993l ). We used a modified Monte Carlo code to simulate the velocity 
distribution expected if the secondary masses were draw n from a distribution described by 



Miller fc Scald (119791 ) as parameterized in Equation 5 of lEggleton et al.l (119891 ). 



Figure [TJ)] shows a 3x3 grid of contour probability plots similar to F igure E] except 
that t he ordinate in each panel is Mi ow , the lower mass limit at which the IMiller fc Scalo 
(119791 ) IMF is truncated. The columns and rows of panels show a range of rj„ and r out , 
labeled in each panel, as in Figure [6] The peak probabilities marginally exceed 60% in four 
panels: middle-left, center, upper-middle, and upper-right. Type Ib/c supernova fractions 
fib/c =0.20 - 0.40 are broadly consistent with the data. However, in all four panels, the 
Opik's Law distribution of separations shown by the solid horizontal line requires a lower 
mass cutoff M\ ow > 1 — 3 M & . The probabilities that the simulated data agree with the 
observations drop to below 40% at Mi ow < 1 M for (3 < 0. The probabilities at low Mi ow 
are somewhat higher in the left column (rj n = 0.02 AU), but we argue previously that 
such small inner limits, coupled with such small (3 are highly improbable physically. Binary 
fractio ns F < 1 would only e xacerbate the problem, requiring even larger lower mass cutoffs 



to the IMiller fc Scalo 

a 



19791 ) IMF. We conclude that the data are grossly inconsistent with 



Miller &: Scald (119791 ) type IMF extending below 2-4 M for any plausible distribution of 



separations. 
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We tested whether a iMiller fc Scald (119791 ) type IMF might be made more consistent 
with the data by the addition of a "twin" component in the same manner described previ- 
ously. Figure [TBI shows a 3x3 grid of contour probability plots similar to Figure [T5] except that 
25% of the stars are now "twins" with q =0.95-1.0 with mean separations < r^ win >= 3.9 AU 
and a log-normal distribution having width <7i ogr = 0.6. The highest contours peak in excess 
of 60%, but these regions lie above lower mass cutoffs of Mi ow =0.8 M e for plausible values 
of r in > 0.063 AU. The Opik's Law distribution of separations would require Mi ow > 0.5 
Mq where it intersects the peak contour rid ge. Even if 25% of th e companions are twins, 
a secondary star mass function similar to the IMiller fc Scald (119791 ) IMF must have a lower 
mass cutoff at relatively large masses in order to produce the velocity distribution in the 
data. The implied type Ib/c supernova fractions are //&/ c >0.35 for the most plausible radial 
bounda ry limits (middle and upper-middle panels), at the upper end of currently measured 
ranges (ILi et al.l 120071 ; iLeaman et al.l 120071 ) . 



We tested the IMiller fc Scald (119791 ) type I MF against the data with a "t win" fraction 
of Tf rac = 0.40, similar to the prescription of iPinsonneault fc Stanekl (120061 ). Figure [T71 
shows a 3x3 grid of contour probability plots similar to Figure [IB] for this series of Monte 
Carlo comparisons. The highest contours exceed 60% over a wide area and approach 70%. 
In nearly all panels, they extend down to cutoff masses of ~ 0.1 Mq. In particular, the 
four upper-right panels having rj n > 0.063 AU and r out > 1000 AU show that the Opik's 
Law distribution of separations i ntersects the ridge of peak contours in the range 0.1 - 0.3 
Mq, i.e., in the range where the IMiller fc Scald (Il979l ) type IMF is flat, or turns over for 
the lowest mass stars. The implied type Ib/c supernova fractions are fib/ c >0.35 for these 



pane l s, at the upper end of, but consistent with, the currently measured values (ILi et al. 



20071 : ILeaman et all 120071 ). 



4.5. Limitations of the Modeling 

In order to compare the observations to the modeled velocity distributions we have made 
a number of simplifying assumptions. Some of these assumptions may account for the result 
that the best models agree with the data at only the 60-70% level. Nevertheless, we contend 
that none of these assumptions affect the major conclusions. 

We assume that all orbits are circular. Although this is unlikely to be the case for the 
majority of systems, it provides a reasonable way to estimate the mean properties of an 
ensemble of systems with a minimum of free parameters, especially in the absence of well- 
sampled velocity curves that would be required to measure eccentricities. It is also a reason- 
able approximation for close, massive binaries where the orbits are likely to have circularized. 
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Early- typ e systems with perio ds less than several days are observed to have eccentricities 
near zero (IGiuricin et al.lll984j ) as are systems with la rge fractional radii, R*/r > 0.24 where 



r is the se mi- major axis and R* the stellar radius (INorth fc Zahnl 120031 ; iPan et al.l 11998 



Zahnlll977l ). Low-mass binary systems have eccentricities near zero for periods shorter than 
~ 11 days (IMeibom fc Mathieull2005l ). Without additional data to provide secure periods 
for much for our sample, we can only say that at least ~ 10 systems appear to be in such 
short period orbits. Furthermore, although the velocity amplitude for a primary star in an 
eccentric system is larger by a factor a/(1 — e 2 ) compared to a circular system of the same 
semi-major axis, the rms velocity dispersion averaged over the orbit is smaller by the same 
factor because the relative radial velocity variations are small over most of the orbital period. 
It follows that the systematic effects of eccentric orbits cancel to zeroth order and can be 
neglected unless the majority of systems are highly ecentric. 

We assume that each primary star has a single companion that dominates the observed 
primary kinematics. In reality, triple and quadruple systems probably exist among the 
sample, but they are likely to be a small fraction of all systems. Among solar type stars the 
fraction of triple and quadruple systems is < 4% (jDuquennoy fc Mayorlll99ll ). The tertiary 
components, if any, are statistically likely to have wide orbital separations and long periods, 
making their observable dynamical influence minimal over few year time baselines and at 
~10 km s _1 velocity precisions. 

We assume that the measured velocity variations are due to orbital dynamics of the 
primary star. Stellar photospheric line profile variations may be present among some of the 
most massive stars, especially the evolved stars in our sample (31 of 114 or about 27% are 
post-main-sequence stars). Line profile variations attributed to atmospher i c pulsations are 
obser ved in >77% of evolved O stars and in some B e stars (jPenrodlll986l ; IVogt fc Penrod 
19831 ) but rarely among dwarfs (IFullerton et al.lll996l ). These phenomena could mimic the 
effects of bona fide orbital velocity variability in systems observed with sufficiently high spec- 
tral resolution to see these effects. Neglecting this possibility would lead to an overestimate 
of the binary fraction and bias the conclusions in favor a larger a and/or smaller (3. Since 
our sample consists predominantly of dwarfs and non-Be stars (there are 3-5 Be stars), the 
impact of line profile variations is minimal and is unlikely to dominate the results. 

Finally, we caution that the binary characteristics within a single OB association may 
not be representative of the binary properties in other massive starforming regions. The 
fraction of binaries, the distribution of mass ratios, the initial orbital separations, and their 
eccentri cities may depend upon the global conditions in the molecular cloud from which they 
formed. I Garcia &: Mermilliodl (120011 ) highlight a possible trend in the binary frequency with 
cluster density and richness. For example, the dense, massive molecular clouds with densities 
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n, 



> 10 7 c m 3 which can produce a super star cluster with M > 10 6 M G) within a diameter 



of a few pc (lO'Connell et al.lll994j ; lElemgreen &: Efremovll 19971 ; iKobulnicky fc Johnsonlll999 



Billett et al.ll2002l ) may produce a different population of binaries than those which generate 



the relatively diffuse OB associations. This would be especially true if binaries result pri- 
marily from gravitational encounters between cloud cores and/or massive protostars early 
in the star formation episode. On the other hand, it might be argued that the production 
of binaries is primarily a local effect driven by physics on the scale of an individual cloud 
core (< 10 AU). If massive stars form primaril y throu gh mergers of intermediate mass stars 
in a dense core as suggested by iBonnell et al.l (119981 ). then perhaps characteristics of mas- 
sive binaries are set by the mass es and separations of subclumps that fail to complete the 
merger process. iKrumhola (120061 ) provides a review of the predictions and consequences of 



competing theories for massive star formation. 



5. Conclusions 

We have analyzed the radial velocity data spanning 6 years on a large sample of OB 
stars from the Cygnus OB2 Association to measure the probable system characteristics 
among massive binaries. We have compared the radial velocity amplitudes to expectations 
from Monte Carlo models based upon several popular formulations for the mass ratios distri- 
butions in massive binaries. These comparisons highlight the allowed parameter space and 
the sensitivity of the models to the key binary parameters: the true binary fraction, F; the 
power-law slope a of the mass ratio distirbution over the range qi ow < q < 1.0; and (3, the 
power law slope of the orbital separations between the inner and outer radial boundaries, 
log(r in ) and \og(r out ). The data improve upon prior studies by including a larger number 
of objects in a common association over more observational epochs. The sample also con- 
sists of a young population with a (mostly) uniform age of 2-3 Myr, thereby minimizing 
evolutionary and selection effects inherent in many samples. 

1. Descriptions of the massive binary mass ratio distribution and orbital separation distri- 
bution derived from spectroscopic and direct imaging studies are biased by selection ef- 
fects and measurement limitations. These produce apparently conflicting prescriptions 
of massive binary parameters. Power-law parameterizations require explicit boundaries 
over which the power law is valid. These boundaries are often insufficiently well speci- 
fied in the published studies to allow meaningful comparisons given that the conclusions 
are sensitive to the choice of r in and r out , as shown in the Figures herein. 

2. Radial velocity amplitude data, such as we model here, are consistent with a range 
of mass ratio distributions and orbital radius distributions. Although the models are 
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degenerate in a and (3, the best fitting power-law slopes are correlated. Under the 
assumption that the distribution of orbital separations is approximately flat (Opik's 
Law = (3 = 0), the data allow only a narrow range of a for given inner and outer radial 
boundaries. 

Power-law parameterizations of f(q) and f(r) for plausible boundary values q min = 
0.02, 0.063 AU < r in < 0.200 AU, r out = 10000 AU, and F = 1.0 are consistent with a 
mass ratio distribution described by —0.6 < a < 0.0 assuming that (3 — 0, i.e., the flat 
Opik's Law distribution of separations. This is consistent with a ~ —0-4 reported for 



adap tive optics imaging studies of Sco OB2 (IKouwenhoven et al.ll2005l ; iKouwenhoven 



20061 ). Best-fitting Monte Carlo Models are, on average, consistent with the data at 
the P a =60-70% level. Agreement between models and data is limited by discrepancies 
between the cumulative distribution functions of the model velocity amplitudes and 
the data both at small velocities Vh ~ 20 km s _1 and at large velocities Vh ~ 80 km s _1 
in roughly equal proportions. This result suggests that small sample size, measurement 
uncertainties, and the possible presence of higher-order systems in the data all act to 
preclude better agreement. 

4. Power-law models with binary fractions F = 0.8 fit the data slightly more poorly than 
models with F = 1.0, and they require larger mean mass ratios with —0.4 < a < 1.0. 

5. Power-law models with binary fractions F < 0.6 fit the data significantly less well, 
having P a < 0.4. They also require extreme values of a, (3, and /j&/ c < 0.20 that are 
inconsistent with observational studies. 



6. The iHogeveenl (119921 ) distribution of mass ratios (a = —2 for systems with q > 0.3 
and a = for systems with q < 0.3 ) for r in = 0.063 AU and r out = 1000 AU can 
reproduce the data at the ~ 60% level, comparable to t he best-fitt i ng po wer-law models 
with a ~ —0.6. This is because the two-component IHogeveenl (119921 ) q distribution 
effectively mimics a power law with slope intermediate between —2 and 0.0. Given the 
ad hoc nature of the this formulation, its additional parameters, and its inability to 
produce a superior match to the data, we find no reason to prefer it over a simpler 
power-law description for f(q). 



7. Secondary star masses drawn from a iMiller fc Scald (jl979j) IMF are inconsistent with 
the data unless the IMF is truncated at low masses below ~2-4 M . IMF truncation 
at such high masses is required to produce the levels of velocity variability seen in the 
data. The implication here is that the physical conditions under which massive stars 
form are conducive to producing relatively close binaries with correlated component 



masses. 
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8. Mass ratio distributions involving both a twin component comprising 25-40% of mul- 



tiple systems and a power-law component with f 



of multiple systems (IPinsonneault fc Stanekl 120061 ) are inconsistent with the data un- 



q) oc q° (i.e., flat) for the remainder 



less the true binary fraction is F < 0.60. Such low binary fractions are shown to be 
improbable. However, mass ratio distributions involving both a twin component com- 
prising ~40% of multiple systems and a power-law component with a = -2 — -3 (i.e., 
Salpeter-like) are consistent with the data. 

9. Secondary star mass distributions comprised of two components, a power-law com- 
ponent and a "twin" component with q > 0.95%, provide good agreement with the 
data at the level of ~ 70%, marginally better than the best-fitting single power-law 
models. If the fraction of twin systems is Tf rac = 25 — 40% and these secondaries have 
mean orbital separations < log[r(AC/)] >= 0.6 ±0.6, then the data are consistent with 
the remainder of the companions being drawn from the canonical Miller-Scalo -type 
field star IMF. This formalism simultaneously predicts a type Ib/c supernova fraction 
°f fib/c =0.30-0.40, consistent with the most recent estimates. This two-component 
description may allow for the production of massive binary phenomena such as the 
large number of dou ble neutron star systems that seem to be required by the data (see 



Kaloaera et al\\2007i . for a review) while simultaneously permitting the production of 



low-mass X-ray binaries at the expected rates. 

10. The best-fitting binary parameters also predict that the fraction of core-collapse super- 
novae that are type Ib/c supernovae is least 30% (assuming single stars do not produce 
type Ib/c supernovae). If the actual Ib/c supernova fraction were 15%, the allowed 
binary parameter space would require binary systems where the two components were 
all nearly the same mass (high a) and with orbital separations well above Opik's law 
(/3 > 1). Such values, although possible, seem unlikely. More likely is that the binary 
observations indirectly support the most recent supernova observations suggesting that 
the Ib/c supernova fraction may be as high as 30-40%. 

Observational work that better defines the distribution of binary separations, mass ra- 
tios, and binary fractions will help to reduce degeneracies in the types of model constraints 
used here. Thereby, such data will allow a more precise specification of initial binary param- 
eters. Optical interferometric observations of companions to massive stars in the radial range 
< 20 AU would help verify that the Opik's Law distribution of separations extends to smaller 
orbital distances. Such data may also begin to directly identify "twin" syste ms with typical 



separ ations ~l-few AU. The continuing Cyg OB2 radial velocity survey (IKiminki et al. 



20071 ) will also measure orbital parameters, modulo inclinations, for a sizable sample of mas- 
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sive binaries with a common origin, permitting tighter constraints on the binary separations 
and companion masses. 
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Fig. 1. — Primary star velocity semi- amplitude versus orbital period for binary systems 
with a 10 M & primary. Solid lines indicate the loci of systems with equivalent mass ratios, 
q = M2/M1—I.O, 0.1, 0.01. Dotted lines indicate loci of systems with common orbital 
separations, r = 0.1, 1.0, 10 AU. The light gray region indicates the portion of parameter 
space sampled in this survey. Namely, the data are sensitive to primary velocity amplitudes 
> 10 km s -1 , periods shorter than ~ 2 x 6 years, and orbital separations < 2 x 10 AU. 
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Fig. 2. — Primary star spectroscopic masses versus observed velocity semi-amplitudes, 14 = 
0-5|14*az — Vmin\, for the 114 sample stars. Typical uncertainties are 5-15 km s _1 . 
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Fig. 3. — The distribution of observed velocity dispersions, V rms (solid line), velocity semi- 
amplitudes, Vh = O.SlKnaa; — V m i n \ (dotted line), and the mean velocity uncertainties (dashed 
line) for the sample. 
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Fig. 4. — Left panel: His togram of the velocity semi-amplitudes for the O stars presented in 
Table 3 and Figure 1 of iGarmany et al.l (Il980l ). Note the differences comp ared to our Fig- 
ure [31 namely a relative lack of stars in the velocity range 20-80 km s _1 in the lGarmany et al. 
( 19801) sample. Midd le panel: Histogram of mass ratios, q, for the O type binaries from 



Garmany et al.l (Il980l ). The histogram is consistent with a flat or rising distribution with 
increasing mass rati o. Right panel: Deriv ed semi-major axis distribution for the O type 
binary systems from I Garmany et al.l (119801 ) assuming an average orbital inclination of 60°. 
The dashed line is a best fit Gaussian curve which peaks near log(r) = —0.8 or ~ 0.15 AU 
with a tail toward larger separations. 
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Fig. 5. — Enclosed regions which denote the parameter space of orbital separation and 
primary mass which produce type Ib/c supernovae for 4 different values of the mass ratio, q. 
Type Ib/c supernovae are preferentially made in systems that undergo a common envelope 
phase, so the parameter space is larger for more extreme mass ratios. Above 32 M Q dot, our 
wind prescription makes all stars (whether in interacting binaries or not) lose enough mass 
to become type Ib/c supernovae. 
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Fig. 6. — Probability contours describing the agreement between the the Monte Carlo sim- 
ulations and the radial velocity data for Cyg 0B2 (solid lines) as a function of a and (3 
in each panel, for a range of minimum orbital radii, r in (columns), and maximum orbital 
radii, r out (rows) with a minimum mass ratio of q = 0.02 and a binary fraction of F — 1.0: 
r in = 0.020 AU (left column), r in = 0.063 AU (middle column), r in = 0.200 AU (right 
column), and r out = 10000 AU (upper row), r out = 1000 AU (middle row), r out = 100 AU 
(lower row). The horizontal line in each panel marks the nominal (3 = (Opik's Law) 
distribution of separations. Contour levels probabilities P a = 10, 20, 30, 40, 50, and 60% 
likelihood that the simulations and the data are drawn from the same parent population. 
The 50% and 60% contours are colored blue and red, respectively, in the electronic edition 
of The Astrophysical Journal. The most probable values exceed 60% likelihood and lie along 
a narrow crescent-shaped ridge line. The dotted contours depict the type Ib/c supernova 
fraction, f Ib / c =0. 10-0.40 in increments o f 0.05, predicted by Monte Carlo binary population 
synthesis models ( jFryer et al.lll998l . Il999l ). 
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Fig. 7. — Histograms showing the distribution of velocity semi-amplitudes, Vh, for the data 
(solid line) and a particular set of Monte Carlo simulations with the indicated values of a and 
j3 (dashed line). Each panel is labeled with the probability, P(Vh), that the two histograms 
are drawn from the same parent population. Both combinations of a, (5 yield probabilities 
of P a ~ 0.60 and lie along the ridge line of peak likelihood from the middle-right panel in 
Figure O Either values of /3 ~ 0.10 paired with a ~ 0.4 or (3 ~ —0.20 paired with a ~ —0.6 
provide similarly good matches to the data. 
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Fig. 10. — Probability contours fro m Monte Carlo s imulations similar to Figure [6] but for 
the mass ratio distribution given by iHogeveenl (Il992l ). 
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Fig. 11. — Probability contours from Monte Carlo simulations similar to Figure [TUl but for a 
binary fraction of F = 0.8. 
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Fig. 12. — Probability contours from Monte Carlo simulations similar to Figure [6] but for a 
two-component mass ratio distribution. This Figure shows a 3x3 grid of contour probability 
plots for T/ mc =0.40 (left column), T/ rac =0.25 (middle column), T/ mc =0.10 (right column), 
and < r Twin >= 15 AU (upper row), < r Twin >= 3.9 AU (middle row), < r Twin >= 1.0 AU 
(lower row). In all cases, the radial limits for the power-law distribution of the non-twin 
component are r in = 0.063 AU and r out = 1000 AU. 
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Fig. 13. — Probability contours from Monte Carlo simulations similar to Figure [12] but for a 
binary fraction F = 0.8. 
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Fig. 14. — Probability contours from Monte Carlo simulations similar to Figure [12] but a 
binary fraction F = 0.6. 
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Fig. 15. — Probability contours from Monte Carlo simulations similar to Figure [6] but for a 
distribution of secondary masses described by the Miller-Scalo IMF with a lower mass cutoff 
shown on the ordinate of each panel. 
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Fig. 16. — Probability contours from Monte Carlo simulations similar to Figure [15] for a 
distribution of secondary masses described by the Miller-Scalo IMF but with a twin fraction 

Tfrac = 0.25. 



-47- 




Fig. 17. — Probability contours from Monte Carlo simulations similar to Figure [15] for a 
distribution of secondary masses described by the Miller-Scalo IMF but with a twin fraction 
T frac = 0.40. 



